/*
 * Normal.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
using System.Text;

using dnA.Math;

namespace dnA.Stats.Distributions {
    ///<summary>This class represents a normal distribution.</summary>
	[Serializable]
    public sealed class Normal:Distribution {

		private static readonly double p_low  = 0.02425;
		private static readonly double p_high = 1.0 - p_low;

		private static readonly double[] ICDF_A = {
													  -3.969683028665376e+01,  2.209460984245205e+02,
													  -2.759285104469687e+02,  1.383577518672690e+02,
													  -3.066479806614716e+01,  2.506628277459239e+00};

		private static readonly double[] ICDF_B = {
													  -5.447609879822406e+01,  1.615858368580409e+02,
													  -1.556989798598866e+02,  6.680131188771972e+01,
													  -1.328068155288572e+01};

		private static readonly double[] ICDF_C = {
													  -7.784894002430293e-03, -3.223964580411365e-01,
													  -2.400758277161838e+00, -2.549732539343734e+00,
													  4.374664141464968e+00,  2.938163982698783e+00};

		private static readonly double[] ICDF_D = {
													  7.784695709041462e-03,  3.224671290700398e-01,
													  2.445134137142996e+00,  3.754408661907416e+00};

        ///<summary>Constructs a normal distribution with a mean of 0 and a variance of 1.</summary>
        public Normal():base(0,1,0,0){}

        ///<summary>Constructs a normal distribution with the given mean and variance.</summary>
        ///<param name="mean">the mean of the distribution.</param>
        ///<param name="variance">the variance of the distribution.</param>
		public Normal(double mean, double variance):base(mean, variance, 0,0){}

        ///<summary>Evaluates the probability density function at X.></summary>
        ///<param name="x">value to evaluate the function at.</param>
        ///<returns>the probability density function at X.</returns>
		public override double EvaluatePDF(double x){ 
			return 1/(StandardDeviation*sqrt2pi)*System.Math.Exp(-.5*(System.Math.Pow((x-Mean)/StandardDeviation,2)));
		}

        ///<summary>Evaluates the cumulative density function at X.></summary>
        ///<param name="x">value to evaluate the function at.</param>
        ///<returns>the cumulative density function function at X.</returns>
		public override double EvaluateCDF(double x){   
			double z = (x - Mean) / StandardDeviation;   
			return 0.5 * (1.0 + SpecialFunctions.Error(z/sqrt2));
		}

		//Port of Sherali Karimov ICDF function at
		//http://home.online.no/~pjacklam/notes/invnorm/impl/karimov/StatUtil.java
        ///<summary>Evaluates the inverse cumulative density function at X.></summary>
        ///<param name="d">value to evaluate the function at.</param>
        ///<returns>the inverse cumulative density function function at X.</returns>
		public override double EvaluateInverseCDF(double d){  
			// Define break-points.
			// variable for result
			double z = 0;

			if(d == 0) {
				z = Double.NegativeInfinity;
			}else if(d == 1) {
				z = Double.PositiveInfinity;
			}else if(Double.IsNaN(d) || d < 0 || d > 1){
				z = Double.NaN;
			}
				// Rational approximation for lower region:
			else if( d < p_low ) {
				double q  = System.Math.Sqrt(-2*System.Math.Log(d));
				z = (((((ICDF_C[0]*q+ICDF_C[1])*q+ICDF_C[2])*q+ICDF_C[3])*q+ICDF_C[4])*q+ICDF_C[5]) / ((((ICDF_D[0]*q+ICDF_D[1])*q+ICDF_D[2])*q+ICDF_D[3])*q+1);
			}

				// Rational approximation for upper region:
			else if( p_high < d ) {
				double q  = System.Math.Sqrt(-2*System.Math.Log(1-d));
				z = -(((((ICDF_C[0]*q+ICDF_C[1])*q+ICDF_C[2])*q+ICDF_C[3])*q+ICDF_C[4])*q+ICDF_C[5]) / ((((ICDF_D[0]*q+ICDF_D[1])*q+ICDF_D[2])*q+ICDF_D[3])*q+1);
			}
				// Rational approximation for central region:
			else {
				double q = d - 0.5D;
				double r = q * q;
				z = (((((ICDF_A[0]*r+ICDF_A[1])*r+ICDF_A[2])*r+ICDF_A[3])*r+ICDF_A[4])*r+ICDF_A[5])*q / (((((ICDF_B[0]*r+ICDF_B[1])*r+ICDF_B[2])*r+ICDF_B[3])*r+ICDF_B[4])*r+1);
			}

			if( d > 0 && d < 1) {
				double e = 0.5 * (1-SpecialFunctions.Error(-z/sqrt2)) - d;
				double u = e * sqrt2pi * System.Math.Exp((z*z)/2.0D);
				z -= (u/(1.0 + z*u/2.0));
			}

			return z*StandardDeviation+Mean;
		}

        ///<summary>Make a deep copy of this <c>Normal</c> distribution.</summary>
        ///<returns>A deep copy of this <c>Normal</c> distribution.</returns>
		public Object Clone(){
			return this.MemberwiseClone();
		}

        ///<summary>Returns the Hashcode for the <c>Normal</c> distribution.</summary>
        ///<returns>The Hashcode representation of <c>Normal</c> distribution.</returns>
		override public int GetHashCode(){
			return (int)(System.Math.Log(Mean)*System.Math.Log(Variance));
		}

        ///<summary>Checks if <c>Normal</c> distribution variable is equal to another object.</summary>
        ///<param name="obj"><c>obj</c> to compare present <c>Normal</c> distributionto.</param>
        ///<returns>Returns true if the variable is the same as the <c>Normal</c> distribution.</returns>
		override public bool Equals( Object obj ){
			Normal value = (obj as Normal); 
			if( value == null ) {
				return false;
			}
			if( this.Mean == value.Mean && this.Variance == value.Variance ) {
				return true;
			} else {
				return false;
			}
		}

        ///<summary>A string representation of this <c>Normal</c> distribution.</summary>
        ///<returns>The string representation of the value of <c>this</c> instance.</returns>
		override public String ToString(){
			StringBuilder sb = new StringBuilder("Mean: ");
			sb.Append(Mean);
			sb.Append(", Variance: " );
			sb.Append(Variance);
			sb.Append(System.Environment.NewLine);
			return sb.ToString();
		}
	}
}
